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Abstract 

P ■ 

We present a description of finite dimensional quantum entanglement, based on a study 
^ ■ of the space of all convex decompositions of a given density matrix. On this space we con- 

CN . struct a system of real polynomial equations describing separable states. We further study 

this system using methods of statistical mechanics. As an example, we apply finally our 
i-Sp ' techniques to Werner states of two qubits and obtain a sufficient criterion for separability. 
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Separability is one of the central issues in quantum information theory (see Horodecki et al. 
[1] for a review) in that in a separable density matrix all correlations are of classical origin and 
no real quantum information processing, as based on the presence of quantum entanglement 
of some kind, is impossible. The solution to the separability problem has been proved to be 
^ ! NP-hard [8] and hence every partial solution constitutes an important achievement. Seminal 
corner stones in that direction have been the Peres-Horodecki criterion (2} [7], and entanglement 
witnessing operators [6]. The first method exploits the fact that positive operators conserve 
the positivity of all separable density matrices, whereas some entangled density operators are 
mapped to non-positive operators. The latter approach uses limits for expectation values of 
>• ! suitably chosen witness operators to distinguish between separable and entangled states. A 
systematic analysis of the so called bound entangled states has been initiated by means of 
unextendible product bases (UPB) [3] which in turn also paved the way towards a formulation 
of the separability problem in terms of roots of complex polynomial equations [4]. As far as we 
know, this route has not been pursued any further and in particular no direct test of separability 
via the convex roof extension of a pure state separability criterion has been probed so far. The 
main obstacles for such an approach have their origin in the complications involved in the 
minimization procedure over all decompositions of the density matrix under consideration. A 
proposal in this direction however has been presented by Osborne [9]. In this work we follow 
this route proposing a similar approach for studying the bipartite separability problem in finite 
dimensional Hilbert space H = Ha <S> 7~£b — C m <8> C n encoding the convex roof minimization 
in a way familiar from statistical-mechanics. 

The paper is organized as follows. After a formal definition of the separability problem and 
a short discussion of pure state separability criteria in the next section, we give a geometrical 
view on the space of p-ensembles and a formulation of the bipartite separability problem in 
terms of a set of nonlinear equations in section [3l A mechanical analogy of these equations 
is drawn in section [4] in terms of a Hamiltonian or cost function on a restricted "phase space" 
and constitutes the basis for the statistical mechanical approach presented in section El After 
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presenting a proof-of-principles calculation for two-qubit Werner states in section [6] we draw 
our conclusions and give a short outlook of the presented formalism. 

2 The bipartite separability problem 

In order to formulate the problem, let us recall the following Definition: 

Definition 1 A state g of a bipartite system AB, described by Ha ® Hb, is called separable if 
there exists a convex decomposition of g composed entirely of product vectors: 

N 

Q=/ J Pi \ x i){ x i\ ® \Vi){ViV \Xi) eH A ,\yi) £7~Lb- (1) 
i=l 

A natural problem arises, known as the separability problem: Given a state g, decide if it is 
separable or not. This problem has been proven to be NP-hard (Gurvits |8J) and (a part of) its 
difficulty lies in the fact that a convex decomposition of a given mixed state g into pure states: 

N 

<? = X>l*«><*il (2) 

8=1 

is highly non-unique (see e.g. Bengtsson and Zyczkowski pH]). Thus, the following Definition 
makes sense: 

Definition 2 Unordered collection {pi, i = 1 . . . N of probabilities and vectors satisfying 

(TJ|) is called a g-ensemble of length N. 

In this work we develop the following approach to the separability problem: we propose 
to search the space of all ^-ensembles of a given state g for product ^-ensembles (^-ensembles 
containing only product vectors), by applying one of the existing necessary and sufficient en- 
tanglement tests to each member of the ensemble. We want the test which has the simplest 
functional form — a polynomial. Such a test is provided by the square of generalized concurrence 
(see e.g. Rungta et al. [llj, Mintert et al. [12], Hulpke [1 3J ) : 

Proposition 1 For any vector G Ha <8> Hb one has that: 

c 2 ^):=M\ A -trn A ^nM(^\) 2 >0 (3) 

and the equality holds if and only if \ip) is product. 

This leads to a set of real polynomial equations describing separable states. The resulting 
system is very complicated due to the fourth order of some equations and a large number 
of variables. Our idea is to study it using methods of classical statistical mechanics. The 
motivation is that such methods have proven to be very efficient not only within classical 
mechanics, but also in many other, distantly related areas (for an application to fundamental 
combinatorial problems see e.g. Kubasiak et al. [T4] and references therein). Hence, we first 
develop a mechanical analogy for our system. Then we define a suitable cost function, or 
"energy", introduce a canonical ensemble, and study the resulting partition function. 
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3 The space of ^-ensembles and separability 



Let us begin with describing the space of all ^-ensembles of a given state g. For convenience 
we pass from normalized ^-ensemble vectors to subnormalized ones: '■= y/pl\^i), such 
that g = YuiLi iV'tKV'il' Let us fix an eigenensemble {|e Q )} of g, where all the vectors \e a ) 
correspond to non-zero eigenvalues Aq, of g, a = 1 . . . r, and r := rank(^) is the rank of g. 
Then, all ^-ensembles are characterized by the well known Theorem by Schrodinger [15] (see 
also PH [L7]) 

Theorem 1 Any g-ensemble of length N > r can be obtained from a subnormalized 

eigenensemble {|e Q )} such that p = Y2 a l e ») (e a | through the following linear transformation: 

r 

\ij>i) ■= ^2z ia \e a ), (4) 

a=l 

where the matrix zi a G C is an N x r block of a unitary N x N matrix, and hence satisfies 

N 

y^^Zj^Zjp = 5 af 3. (5) 
1=1 

Theorem [Q gives us the characterization of all possible ^-ensembles in terms of N x r 
matrices z, satisfying the condition j5j). Geometrically, this condition defines the so called 

Stiefel manifold 

V N , r :=U(N)/U(N-r). (6) 

It forms a principal fiber bundle over the Grassmann manifold Gn,t (the set of r-dimensional 
subspaces of C N ) with a fiber diffeomorphic to U(r) (we refer to Kobayashi and Nomizu 
Vol. 1 [18] for the definition and basic properties of fiber bundles and to Spivak Vol. 5 [T9] for 
more information on the Stiefel and Grassmann manifolds). 

However, note that there is some additional symmetry: from Eq. ([2j) we see that the order of 
vectors in a ^-ensemble does not matter, and thus two JVxr matrices z, z' satisfying Eq. (JSJ) and 
differing only by a permutation of their rows define the same ^-ensemble. To fix this freedom, 
observe that a z-matrix satisfying Eq. ([S]) has necessarily rank r, and hence we may consider 
only those matrices z, for which the first r rows are linearly independent. The set of such z J s 
constitutes a simply connected open subset of Vn, t (which is nevertheless dense in Vjv, r ) and 

U(r) 

over such a neighborhood the bundle Vn, t ► Gn, t is trivial by construction. This allows us 

to formally write down an explicit solution of the constraints (jSJ) 

z = GS^yU, (7) 

where U G U(r), l r is the rxr unit matrix, v is an arbitrary, complex (N-r)xr matrix, and GS 
denotes the Gram-Schmidt orthonormalization [20] applied to the columns. There are no more 
symmetries, since we have defined in Definition [2] ^-ensembles using vectors \ipi) rather than 
more physical projectors (ipi\ , as the latter are harder to work with. In case of ^-ensembles 
defined through projectors, there would be an additional symmetry of multiplying each row of 
z by a (different) phase. Comparing Eq. (jZJ) and Eq. ((H), one sees that an arbitrary ^-ensemble 
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of length N is obtained from the fixed eigenensemble by i) applying a unitary rotation to |e Q )'s 
and ii) subsequent increasing of the length of the ensemble along the Grassmannian G^r- 

So far we have characterized ^-ensembles of a fixed length N. It seems that in the search 
for product ensemble we would have to consider all possible lengths N > r. However, from 
Caratheodory's Theorem (see e.g. Kelly and Weiss [21]) we know that a separable state can be 
decomposed into at most iV = m 2 n 2 linear independent (in M m n _1 ) product states. Hence, it 
is enough to consider only ^-ensembles of the length N = m 2 n 2 . (there is a natural inclusion 
of space of shorter ensembles in the space of longer ones). 

Let us now examine the entanglement test given by Proposition Q3 First, we quote some 
well known facts regarding the geometry of pure product states (see e.g. Bengtsson and 
Zyczkowski [10]). Note that the polynomial c 2 (ip), defined in Eq. (|3j) , is in fact a sum of 
modulus squared of quadratic, complex-analytical polynomials in \ip): 

c 2 w = \ E i <cr ® cr\r B ® r B ')\\ ® 

a,b=l 

where {\Ca A ')}a=i,...,d 1 , {\( BB ')}b=i,...,d 2 are orthonormal bases of the skew-symmetric spaces 
U A A H a > = C m AC™ and H B A H B > = C n A C\ respectively. Thus, c 2 {ip) = 0, and hence \ip) 
is product, if and only if 

(Ca® (b\ip® V) = for all a, b. (9) 

It is worth noticing that this is just the condition for the matrix of components of to have 
rank one. Geometrically, the system of homogeneous equations ((9|), or equivalently the single 
equation c 2 (ip) = 0, describes the image of the so called Segre embedding CP m x CP" c — > CP mn 
given by ([x], [y]) h- ► [x <S> y] ■ As we can see from Eqs. j9j), this image, i.e. the set of product 
vectors, is a complex-analytical manifold — as an intersection of complex quadrics — in contrast 
to the Stiefel manifolds Vn,v, which are real. 

Since for alH = 1, . . . N polynomials c 2 (ipi) are non-negative and equal to zero if and only 
if \ipi) is product, we can sum them up for a given ^-ensemble, and thus obtain a collective 
separability test for the whole ^-ensemble, given by a single polynomial function. Combining 
this with the parametrization (j4j) and the constraint (jSJ), we obtain the following description of 
separable states: 

Proposition 2 A states g of rank r on C m ®C™ is separable if and only if the following system 
of equations possesses a solution 

9 9 9 9 

m n m n r 

E e( z ) ■= E = E E z-z^Kp^ z i^ = °» ( 10 ) 

i=l i=l a,...,u=l 

9 9 

m n 

Ca.p{z) '■= ZiaZi/3 — 8 a p = 0, (11) 



i=l 



where 



K^ v ■= \(e a ®ep\Yl m ®Yl ri e^®e v ) (12) 

and n m , n n are the projectors from C m ® C m , C n <g> C n onto the skew-symmetric subspaces 
C m A C m , C n A C n respectively. 
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We note that the pure state entanglement measure we use is the square of the generalized 
concurrence c(g) (cf. Rungta et al. [TT], Mintert et al. [32]), however, instead of the convex 
roof construction c(g) := inf YliPi c (^i) = i^X^ ^) (where are normalized vectors), we 
analyze 

E 8 (z) (13) 

i 

as a "quantifier" of entanglement. We remind the reader that no caveats are introduced by 
this, since we are only interested in the detection of zero entanglement rather than in the full 
construction of an entanglement monotone. 

Note that the Eqs. ( jTOl ) . fTTTT l are invariant with respect to local unitary transformations, 
since when g is separable also Ua <£> UbqU\ ® U B is, for arbitrary Ua £ U(m), Ub £ U{n). The 
latter transformation can be viewed either as a local change of basis (passive view) or as an active 
rotation (active view). Indeed, from Eq. ([3]) one immediately sees that (?{U a ® Ubi/j) = c 2 (ip). 
Thus, the function E g and all quantities derived from it are constant on the whole unitary class 
of g, i.e. on [ e ] ■= {U A ® U B gU A ® U B ; U A £ U(m), U B £ U(n)}. In what follows, we refer 
with g to its local unitary class [g]. 

We give a brief comparison to a previous analysis carried out by Wu et al. in Ref. |22j, 
also leading to a different set of polynomial equations. These authors have used a higher 
order polynomial test for separability: let a a '■= tr^ B |^)(^|, then \ip) is product if and only 
if det(<T4 — 1) = 0. The relation to Eq. ([3l) is established by observing that det(o4 — 1) = 
J2T=o(~ l) fec fc( cr A), where cVs form a basis of [/(m)-invariant polynomials (see e.g. Ref. [23]). 
Particularly, 2c2(o"a) = (trcr^) — trcr^, which is precisely the generalized concurrence squared 
(cf. Eq. ([3])). For testing separability, it is sufficient to consider only c 2 . 



4 Mechanical analogy 

Equations (TTOl) and (TTT1) form a system of real (after taking real and imaginary parts) polynomial 
equations. Let us denote by V e the set of its solutions for a given g. Then the separability 
problem is equivalent to the question whether V e is empty or not. In principle there is a general 
solution to such problem, provided by the so called Real Nullstellensatz (see e.g. Bochnak et 
al. [H] )• It says that = if and only if the ideal generated by the polynomials E g , 
{ReC Q/ 3, ImCo^}, and by all (real) sum-of-squares (SOS) polynomials^ contains the constant 
— 1. Equivalently, V g = if and only if there exist a SOS polynomial s = ^ n (^n) 2 , a real 
polynomial t, and (complex) polynomials u a p such that: 

-1 = s(z) + E e {z)t{z) + Re [Caf}(z) u a p(z)] . (14) 

<x,P 

However, finding such a certificate is computationally very difficult and inefficient, due to the 
fact that the degrees of polynomials s, t, and u a /3 are a priori unbounded (see also Refs. [91125]). 

Here we develop a different approach based on a statistical analysis of a classical-mechanical 
analogy. Namely, we treat z ia as a collection of complex row vectors Zj £ C r = R 2r , i — 1 . . . N 
and treat each row Zj as a complex phase-space coordinate of a fictitious particle moving in 

1 Interestingly, SOS polynomials also appear in a solution to the classicality problem of states of a single 
mechanical system — they are enough to detect a very broad family of states through generalized squeezing 
conditions (Korbicz et al. |26|). 
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r-dimensional space. Then the whole matrix z ia becomes a phase-space coordinate of a system 
of N such particles in their composite phase-space T : = R 2r x • • • x M 2r = C Nr . Now, let 
E e (zx, . . . , Zjy) and C a p(zi, . . . ,z N ) be defined by Eqs. (fTOl) and f fTTI ). We emphasize that E g 
depends on the separability class of the analyzed state g through the fixed eigenensemble {|e Q )}. 
From the property ([3]) it follows that 



We will think of E e as a cost function or Hamiltonian (it is extensive in the number of fictitious 
particles N), of our fictitious mechanical system. Then, we can treat C a p as at the primary 
constraints imposed on the a priori independent phase-space coordinates (zi, . . . , zjy). We 
note that even if the mechanical system corresponds to free particles (if E g was diagonal) the 
resulting model is nevertheless interacting due to the forces of inertia induced by the non-linear 
constraints. 

The corner stones of the mechanical interpretation of the separability problem (fTOl) . (fTTi l 
can be summarized as follows: the ^-ensembles of density matrices with a fixed rank r form 
the Stiefel manifold Vjv >r , which we can be viewed at as a constraint surface in the phase-space 
T. Each state g defines the non-negative cost operator E e a ^ (TT2]) which uniquely defines the 
cost function E e on T, which probes the separability of the ensembles. The cost function E e 
assumes the value zero on the constraint surface Vn, t (which is then its global minimum) if and 
only if g is separable. 

5 Statistical-mechanical approach 

Although in principle one could tempt to solve the constraints explicitly by Eq. (£3), the resulting 
parametrization of the constrained manifold is rather hard to work with due to the iterative 
nature of the Gram-Schmidt orthonormalization. We circumvent the complications with an 
explicit incorporation of the constraints by using a standard method of implicit treatment of 
constrained systems due to Dirac [27]. It is based on the introduction of Lagrange multipliers. 
To this end we define the full Hamiltonian of the systems as 

Hf u u(zi . . . Z N ) := E e + ^^LOapCap, (16) 



where u a p are the Lagrange multipliers. Note that the constraints written in the matrix C a p 
are not all independent: it is in fact a hermitian matrix and we need to employ one Lagrange 
multiplier for each independent constraint only. On the other hand we have considerable 
freedom for choosing the spurious Lagrange multipliers in the Lagrange matrix uj. We choose 
uj to be hermitian. Then, Hf u u is hermitian and has only real eigenvalues. Moreover, in order 
to take into account all independent constraints, we require that detcu 7^ 0. The constraints 
C a p = are then realized on average by setting to zero the variation of Hf u u with respect to 
Map : dHf u u/du a p = 0. 

The number of fictitious particles N will in general be notably large — in dimension 2 ® 4 
for example we have iV > 64. Thus, the direct analytical study of our fictitious mechanical 
system seems rather hopeless and we proceed further using methods of statistical mechanics 
and numerical simulations. The most natural framework would be microcanonical ensemble, 
however it is also difficult to work with. Hence, we will introduce a canonical ensemble, keeping 
in mind that this is just a technical tool, so, for example, the inverse temperature (3 plays only 
a role of a parameter here, without any physical meaning. 



E g (z u . . . ,Zjv) > for any (z x , . . . 



z N ) e C 



■ Nr 



(15) 
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We proceed to define the canonical partition function Z for our system. The most natural 
definition is perhaps the following 

J i,H a<f3 

r JV n , N r 

= / n n s y^.zfoZjp - <w ex p^ - ^ ^2 z-z^e^ 

i,fi c*</3 i=l i=l a,...,u=l 



(17) 



where the integration is explicitly restricted to the constraint surface V^ r (cf. Eq. ([£])) given 
by C a/ 3 = 0. The intuition behind such an approach is the following. We can formally introduce 
constraint "state density" function 



p(e) := ! l[d 2 z lfl l[5[C a p(z 1 ...z N )]5(e-E g (z)). 
Since E g (z) > 0, p(e) is non-zero only for e > 0. Then: 



(18) 



Z{(3- g) 



dep(e)e 



-/3e 



(19) 



Let us assume that the state in question is entangled. Then E e (z) is strictly positive, so there 
exists a constant a such that E g (z) > a > 0. The average "energy" is then separated from zero: 



({E e )) 



Z((3 



| poo j roo 

%F)L iep{ ' )ee '"-m7)L d<p(£,oe 



-0e 



a. 



(20) 



Now let q be separable. Then, by Proposition [2] E e (z) has zeros on the constraint surface 
C a/ 3 = with each zero corresponding to a separable ^-ensemble. Since such ensembles are 
"rare", we expect that the state density p(e) — > with e — >• 0. Let us assume for a moment that 
the leading term in the actual dependence of p(e) was given by a power law 



p(e)=Ae s , A,5>0. 
Then we obtain the well established result Z(j3; g) 
T(5 + 2) 1 5 + 1 



T T(5 + 1) and 



({E e )) 



r(8+l) (3 (3 



(21) 



(22) 



Thus, we put forward the following conjecture: 

Conjecture. For the "state density" function p(e), defined in Eq. ( figjj . it holds: i) the 
mean energy (defined in Eq. fSW ) ((E e )) = a > if and only if g is entangled; ii) the mean 
energy ((E Q )) scales as 1/(3 if and only if g is separable. 

We anticipate that indeed we observe such a behavior in a simple case of 2 <g> 2 Werner 
states [28]. Note that in general the exponent 6 will depend on the state 5 = 6(g). 

The partition function defined in Eq. (Tl~7l) is difficult to work with analytically (however 
one can still investigate it numerically, e.g. using Monte Carlo methods), so we use a different 
object — the partition function for the full Hamiltonian ( 1T61) . We first rescale the variables: 
Zia i— > Zia/yN and then define: 



Z{(3,u; g) :-- 



IF 



i, 1 1 



z ifl exp 



[E e {zi ...z N ) + Nj2( z i\^i) - Nhvcu 



N 2 



(23) 
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where ( • | • ) denotes the standard scalar product in C r . Performing further rescaling: 



Z(/3,uj; g) becomes (after dropping the tildes): 



Z(p,u; g) 



/ JJ d 2 z ifl exp J - (3E e (zi . . . zjy) - ^(ziluzi) + Ntxu 

i,fj, i 



Y[ d 2 z itl exp - ^2 Zia Zi P E a 



— ^^(zj|ciJZj) + iVtra; 



Now we are able to reproduce the (rescaled) constraints (fTTjl only on average: 




a/3 



\ogZ(P,u; q) = ((N5 af3 -^ZiqZjp)) 



(24) 



(25) 



(26) 



where the average ((•)) is taken with respect to the probability density defined through Eq. j25 

1 



P e (zt . ..z N ;/3,u) 



exp 



Z(j3,u\ q) 

Thus, requiring that d/duj a p\ogZ(/3,u; g) — amounts to: 
N5 a/3 = {(^z—Zip)). 



(3E e {zi . . . z N ) - y ^z^uzj) + Ntxu 



(27) 



(28) 



Following the standard treatment of constrained systems, the equations ( 1281) are treated as 
conditions imposed on a priori arbitrary (apart form being hermitian and non-singular) matrix 
of Lagrange multipliers u>. We note that the above approach based on Hf uU is nothing else but 
a (formal) evaluation of the integral (TT7l) through the saddle point method with N — > oo. 

A significant simplification of the partition function (l25l ) comes from the form of our Hamil- 
tonian E e — from Eq. ( jTOl) it follows that E e [z\ . . . z^) = ^2 i Ei g (z i ), where E e is just the 
function E g with N — 1. The situation is more subtle with the constraints ( |28l ). For the pur- 
pose of this work we will assume that the contribution to the sum from each fictitious particle 
is equal, i.e. ((2^2^) ) = <5 Q/ 3 for every i. In general, such "equipartition" of course does not have 
to hold and it is an additional restriction on the Lagrange multipliers. By such an assumption 
we however achieve a factorization of the partition function: 



Z(/3,u>;q) = [Z 1 (J3,u;q)] n , 
where Z\ is a one-particle partition function: 



(29) 



Y\ d% exp [ - PE g (z) - (z\uz) + tru 



J\ d% exp |^ - (3 ^2 z » z P E lp»v z » z v ~ ( z k z ) + trw 

fj.=l a,...,v 



(30) 
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From now on we will consider Z\ only. The constraint equations (l28l) are then replaced by 
a one-particle version: 







duj, 



a/3 



logZi(/3,o;; g) = 5 a p - ((z a zp)) = 



(31) 



in accordance with our extra assumption made above. The average in Eq. f|3TI) is taken with 
respect to the probability distribution: 



Zi(P,u; g) 



exp 



(3E\ g (z) — (z\uz) + tru; 



(32) 



In particular, Eq. (IBTl ) implies that ((\z a \ 2 )) = 1. 

To understand the meaning of Eq. ( T3TT ) . let us assume that u = uio{/3] g) is its solution. Then 
Eq. (I3~TT) implies that a family of vectors {\ip{z)) := ^2 a z a \e a ) ; z G C r } forms a continuous 
^•-ensemble with respect to the probability distribution (l32l . i.e.: 



d*zP lfl (z;^a;o)|V(z)X^(z)| = e 



(33) 



irrespectively of (3. Since £i e (z) = c 2 (ip{z)) (cf. Eqs. © and (H)) is the concurrence squared 
of each \ip{z)) , the average "energy" is just the ensemble average of the concurrence squared: 



/r\ 
d 2j zP le [z;/^ (/3; Q)] Ey e {z) = "^log^i 



w=u>o(P;q) 



(34) 



Due to the property (jHJ) one can formally simplify the integral ( 1301) using the Hubbard- 
Stratonovitch trick. Indeed, Eq. (l30l) can be rewritten as: 



Zi(P,u; g) 



/r , di,di 

Ai=l ^ a,6=l a,/3 



(z|a>z) + tra> >, 



where: 

h%{0) ■= (Ca ® Cfc|e« ® e^) 
and we have rescaled (3 by 1/4. Next, we use the Hubbard-Stratonovitch substitution: 

]2„ / |„|2 



exp(-/3|y| 



d z s /| s 

— exp(-— + isy + isy 



(35) 



(36) 



(37) 



to obtain (after a formal interchange of the integrations) : 

Zi{p,u\ g) = 



d\,d2 .o 

j-j- d-.s,,/, 



exp 



a,b=l r 

^2 [ ~ z^u aj3 zp + i ^ Sg~bhu (g)z a Zf3 + i ^ s afc h a ^(g)z^z^ 



xexp<j > I - 

a,/8 



(38) 



The above integral is finite if and only if uo > (as we said earlier we assume uo to be non- 
singular in order not to loose any of the constraints, hence the strong inequality here). This 
puts no restriction on the amount of independent parameters in u and from now on we will 
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assume this condition to hold. Performing the Gaussian integration in the 2r variables z, z 
finally yields: 



Z 1 ((3,u; 6 ) = 7T r J j|±M eX p(-I^ 

a, 6=1 a, ft 



\s ab \ + tTLV 



^detM e (s,uj) 



(39) 



where 2r x 2r matrix M Q (s,u) is defined as follows: 



M 8 {s,uj) :-- 



-2iE a , 6 ^ ha6 (^) S7 



(40) 



(we used the fact that u = uj t ) and h ab (g) denotes the r x r matrix whose elements are h^^g). 



6 Calculation for Werner states 



In this Section we apply the developed statistical method to study Werner states of a 2 <g> 2 
dimensional system. They are defined as follows: 



W(p) := (l-p)|*_)<*_| + |l 2 ®l 2 , 



where: 



l*±> == ^(I01> ± |10», |<&±> := ^(|00) ± |11» 



(41) 



(42) 



are the Bell basis states and {|0), |1)} is the standard basis of C 2 . The states W(p) have positive 
partial transpose, and hence are separable (Peres and Horodecki et al. [2j), for p > 2/3. As 
the fixed eigenensemble {|e Q )} of W(p) we take: 



\ei) 



l-^l*->, |e 2 ):=^i|^ + ), 



M := f |* + >. 



(43) 
(44) 



We proceed to calculate the one-particle partition function Zi\J3, u\ W(p)) = Z-i(P,uj;p). In 
what follows we assume p > 0, for p = corresponds to a pure state. According to the general 
formula (l39l) . we have to find the matrices h ab (W(p)) and M w ^(s,u), defined in Eqs. (|36l) and 
( 1401 . Since in the case of C 2 <g> C 2 the skew-symmetric subspace C 2 A C 2 is one-dimensional — it 



is spanned by a single vector \() = l/y2(|01) — |10)) in each copy AA' and BB' — there is 
only one matrix h ab (W(p)) = h(p) and only one Hubbard-Stratonovich parameter s a b = s. 
Calculation of h(p) and M W ( p ){s ) uj) = M p (s,uj) yields: 



Hp) = o 



M p (s,u) 



4-3p 



p 

p 

p 

uj — 2ish(p) 
— 2ish(p) uj 



(45) 
(46) 
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so that: 



1 2 

Ei(z;p) = — 1(4 - 3p)z 1 + pz\ + pz\ + , 



(47) 



and: 



Zx(/3 , u; p) — I d 2 zi . . . d 2 z^exp — (3 I (4 — 3p)z'( + pz\ + pz 2 + 1 2 — (z\uz) + tvu (48) 



(we have absorbed the factor 1/64 into the definition of the parameter (3). 

Next, we calculate detM p (s, uj) for p ^ 0. We first perform a transformation: 



M p ^ :-- 



hip)- 1 / 2 
h(p)-V2 



M„ 



hip)- 1 / 2 
h(p)- 1 / 2 



uj' -2is 
-2is H' 



where: 



uj' := hip)- 1/2 uhip)- 1/2 . 
Then we multiply Eq. (l49~I) on the left by 



(49) 



(50) 



1 

2is uj' 



to obtain: 



1 

2is uj' 



M' 

p 



UJ 



-2is 



4|s| 2 + uj' uj' 
and after taking the determinants of both sides: 
detM p is,u) = deth(p) 2 det(4|s| 2 + u'u7). 
We then substitute Eq. (l52l into Eq. (l39l and finally obtain (with x := 4|s| 2 ): 



(51) 



(52) 



Zi(/3,u;p) 



7T 



4/5 deth(p) 



3 tr[a/h(p)] 



dxe *0 



det (a; + cj' ct/) 



(53) 



where uj' is defined through Eq. (IBTil) . The above integral is well defined, since det ( x + uj' uj' ) = 
det(x + Vu? uj'\fuj') and \/ r u? uj' y/u? is strictly positive, as we have assumed that detuj ^ 0. 
We can explicitly calculate the derivative dlogZi{/3,uj;p)/duj'. For a generic uj' it takes the 
following form: 



<91ogZi(/3,u;;p) 



duo' 



<91ogZi(/3,u;;p) 



VW) 



dye 



-i —1 oo 



4/J 



dxe 4 f 



det (y + uj' uj') J ^ J det (x + a;' uj') 



h(p) — [x + uj' uj') 1 uj' 



(54) 



The special case of Eqs. (15311 . (15411 for Bell-diagonal states is straightforward — it is enough to 
replace matrix h(p) from Eq. (I45ll with the diagonal matrix 4diag(l —p\—pi — PziPiiP2-,Pz)- 
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Figure 1: The plot of min | \d\ogZ 1 (P, uj\ p)/duj'\\ HS for Werner states as a function of probability 
p for (3 = 10. 



7 Numerical results 



Further studies of the integral (1531) were performed using numerical methods. According to 
Eq. ( 13T1) one has to search for a saddle point of logZi(/3, uj;p) with respect to uj (or equivalently 
with respect to uj'; cf. Eq. (1501). The search was performed by flood-minimizing the Hilbert- 
Schmidt norm of d\ogZi(j3 , uj; p) / duj' a/3 for a range of parameters (3 = 10, 100, .... For simplicity 
we assumed a specific form of uj': 



UJ 



7 

A 

A 

A 



(55) 



and minimized the derivative (given by formula similar to to Eq. ( 1541 ). but taking into account 
the specific symmetry of (|55l) ) with respect to the parameters 7, A > 0. We payed attention that 
the obtained minima are not on the border of the region uj' > (or equivalently uj > 0). The 
specific choice (l55l) of uj' was motivated by the form of the cost function (1471) . We also obtained 
some numerical evidence that in the generic case the minima of \\d\ogZi(j3,uj;p)/duj'\\Hs were 
attained for matrices uj' very close to ( |55l) . The results of the simulations for /? = 10 are 
presented in Fig. Q] (the results for higher values of (3 did not differ from those for (3 = 10). We 
see that for p > 0.89 the constraints (IBTT) can be satisfied. We shall call the interval where it 
happens "equipartition region". 

Next, the dependence of the average entanglement ((E lw ^)) Q (f3) (cf. Eq. (1341 )) of the 
continuous ensemble (1331) on j3 within the equipartition region was examined (recall that outside 
this region the one-particle constraints (|33l) are no longer satisfied). Fig. [2] shows a sample plot 
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Figure 2: The plot of ({Eiw(p)))o{ft) for p = 0.90 on a double log scale. 



for p = 0.9. One sees that the average "energy" indeed scales like 1//3, just like predicted by the 
Ansatz l j2~T| ) and Eq. ([221 . The estimated exponent 5 at this value of p is 5 ~ 1.75. We have 
also checked that in the limiting case [3 — > oo the equipartition region is not altered. Hence, 
our procedure seems to detect separability of the Werner states ( 1411 at least for p > 0.89 and 
thus can serve only as a sufficient condition for separability. We did not check the behavior of 
((Eiw(p)))o(^) outside the equipartition region p < 0.89. 

8 Further questions and concluding remarks 

The statistical mechanical approach to the separability problem as presented here differs from 
the more traditional techniques in that we studied the space of convex decompositions of a 
given state, rather than the convex set of all states. The resulting polynomial equations are 
real due to the constraint (fTll and this real structure makes the analysis more complicated 
than it would be in a complex case. Hence, we applied statistical-mechanical methods to study 
possible zeros of this system. As an example we studied 2 <g> 2 Werner states (l4Ti ). However, 
the numerical difficulty already at this simple example was quite high and we have applied 
several simplifications. Nevertheless, the numerical results suggest that at least for separable 
states in a vicinity of the identity, the partition function and the average "energy", related to 
the ensemble entanglement (cf. Eq.(fT3l). show some qualitative change in their behavior. 

There are obviously some important questions left. First of all, we postulated rather than 
derived the power-law state density behavior (1221 ) for entangled states. It would be an inter- 
esting, albeit difficult, task to try to analytically derive this law. Or at least to find some 
arguments in its favor. 

Another thing is that in passing from the full iV-particle constrains (l28l) to the one-particle 
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one (l3~TT) we have tacitly assumed a sort of "equipartition" of the constraints, i.e. that the 
constraints are divided equally among the particles. But it actually does not have to be like 
that. In particular, the shape of the curve in Fig. [Q tells us that below p = 0.89 the constrains 
are not "equiparted". Thus, in principle one should work with the full A-particle partition 
function (l25j) and seek regions were full constraints ( )28"]) can be satisfied. Then the scaling of 
the average "energy" with /3 within that regions will be able to discriminate between separability 
and entanglement. 

As a side remark, we note that quite surprisingly, the value p = 0.89 appears in Braunstein 
et al. [29] separability criterion, based on an estimation of the size of a ball of separable 
states around the normalized identity (see also Bengtsson and Zyczkowski [10] and the refer- 
ences therein). It will be worth analyzing this curious coincidence in order to gain a deeper 
understanding of the strengths and weak points of the presented approach. 

Finally, let us mention that in principle one can try to directly numerically calculate integral 
lfl7|) using Monte Carlo method. The points of Vn, t can be generated either using Eq. (JZj) or, 
what seems more feasible, directly from definition ([6]). The latter method amounts to generating 
random unitary matrices from U (N) and discarding (N — r) of their columns (for the methods 
of random generation of unitary ensembles see e.g. Pozniak et al. [30]). However, we have 
not performed such simulations. 

We gratefully acknowledge discussions with H.-U. Everts, P. Horodecki, G. Palacios, and 
R. Wimmer. We would like to thank Deutsche Forschungsgemeinschaft (SFB 407, SPP 1078, 
GK 282, 436 POL), the European Graduate College 665, EU IP project SCALA, ESF PESC 
Program QUDEDIS, Spanish MEC Program Consolider Ingenio 2010, and Trup Cualitat Gen- 
eralitat de Catalunya for the financial support. 
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